Si dos variables \(X\) y \(Y\) tienen covarianza \(S = \begin{pmatrix} a && b \\ c && d \end{pmatrix}\), muestre que cuando \(c \neq 0\), la primera componente principal está dada por \[ \sqrt{\frac{c^2}{c^2+(V_1-a)^2}}X+\frac{c}{|c|}\sqrt{\frac{(V_1-a)^2}{c^2+(V_1-a)^2}}Y \] donde \(V_1\) es la varianza explicada por la primera componente principal.
Considere los datos en el archivo T8-5.DAT correspondientes a un tramo censal. Suponga que las observaciones en la variable \(X5\) = valor de la mediana de hogares fue registrada en unidades de diez miles más que de cientos de miles de dólares (es decir, multiplique todos los datos listados en la sexta columna por 10).
En cada caso, compare las estimaciones con los datos en diez miles y cientos de miles (son dos matrices de covarianzas) para las componentes pincipales, trate de obtener una interpretación y de explicar el efecto de cambiar de escala.
col_names_censo <- c('poblacion', 'p_profesional', 'p_empleados_m16',
'p_empleados_gob', 'mediana_hogares')
censo <- read_table('./data/T85.DAT', col_names=col_names_censo)
rm(col_names_censo)
censo_escala2 <- censo %>%
mutate(mediana_hogares = 10*mediana_hogares)
pc_orig <- prcomp(censo, scale=F) %>%
print()
## Standard deviations (1, .., p=5):
## [1] 2.6326932 1.3360929 0.6242194 0.4790918 0.1189747
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4
## poblacion -0.78120807 0.07087183 -0.003656607 0.54171007
## p_profesional -0.30564856 0.76387277 0.161817438 -0.54479937
## p_empleados_m16 -0.33444840 -0.08290788 -0.014841008 0.05101636
## p_empleados_gob -0.42600795 -0.57945799 -0.220453468 -0.63601254
## mediana_hogares 0.05435431 0.26235528 -0.961759720 0.05127599
## PC5
## poblacion 0.302039670
## p_profesional 0.009279632
## p_empleados_m16 -0.937255367
## p_empleados_gob 0.172145212
## mediana_hogares -0.024583093
summary(pc_orig)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 2.6327 1.3361 0.62422 0.47909 0.11897
## Proportion of Variance 0.7413 0.1909 0.04168 0.02455 0.00151
## Cumulative Proportion 0.7413 0.9323 0.97394 0.99849 1.00000
pc_nuevo <- prcomp(censo_escala2, scale=F) %>%
print()
## Standard deviations (1, .., p=5):
## [1] 7.1392827 2.5787380 1.1916363 0.4793601 0.1190091
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4
## poblacion 0.05765977 -0.78220765 -0.02307185 0.541282426
## p_profesional -0.03281610 -0.35032828 -0.76422159 -0.540439691
## p_empleados_m16 0.03467366 -0.32717525 0.10109097 0.051177198
## p_empleados_gob 0.07557524 -0.39085238 0.63207658 -0.642117054
## mediana_hogares -0.99432619 -0.07491367 0.07545108 0.002204195
## PC5
## poblacion -0.302171422
## p_profesional -0.009137870
## p_empleados_m16 0.937504990
## p_empleados_gob -0.172301138
## mediana_hogares 0.002375237
summary(pc_nuevo)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 7.1393 2.5787 1.19164 0.47936 0.11901
## Proportion of Variance 0.8598 0.1122 0.02395 0.00388 0.00024
## Cumulative Proportion 0.8598 0.9719 0.99589 0.99976 1.00000
En la escala original, la primera componente principal explica el 74% de la varianza. Le asigna el coeficiente de más peso a poblacion, pero salvo por mediana_hogares, las otras variables tienen pesos moderados y similares.
Al hacer el cambio de escala en mediana_hogares, sin embargo, la primera componente se vuelve casi mediana_hogares por si sola, pues los demás coeficientes son casi cero. Esto puede explicarse porque componentes principales está pensado para explicar la mayor varianza posible, y en esta nueva escala, como \(\textrm{Var}(10X_6)=100\textrm{Var}(X_6)\), la varianza aumenta significativamente, y también lo hace en comparación con las demás variables, como podemos ver usando
apply(as.matrix(censo), 2, sd)
## poblacion p_profesional p_empleados_m16 p_empleados_gob
## 2.0754652 1.3294632 0.8948008 1.4033798
## mediana_hogares
## 0.7101973
apply(as.matrix(censo_escala2), 2, sd)
## poblacion p_profesional p_empleados_m16 p_empleados_gob
## 2.0754652 1.3294632 0.8948008 1.4033798
## mediana_hogares
## 7.1019731
Consideren la matriz de correlaciones dada abajo. Los datos originales corresponden a las mediciones de 8 variables de química sanguínea de 72 pacientes en un estudio clínico. (Jolliffe, 2002). La matriz de correlaciones de las variables rblood, plate, wblood, neut, lymph, bilir, sodium y potass, en ese orden, está dada en \(S\), y la desviación estándar de cada variable en \(\mathbf{\sigma}\).
Aplique componentes principales a la matriz de covarianzas y a la matriz de correlaciones y explique las diferencias. Basado en la observación anterior, ¿sobre qué debería hacerse el análisis?
R <- read_csv('./data/5-S.csv')
sigma <- c(0.371, 41.253, 1.935, 0.077, 0.071, 4.037, 2.732, 0.297)
names(sigma) <- names(R)
S <- t((t(as.matrix(R) * sigma)*sigma))
pc_S <- princomp(covmat=as.matrix(S), cor=F)
pc_R <- princomp(covmat=as.matrix(R), cor=T)
pc_S$loadings
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## rblood 0.943 0.329
## plate -0.999
## wblood -0.192 -0.981
## neut 0.758 0.650
## lymph -0.649 0.760
## bilir 0.961 0.195 -0.191
## sodium 0.193 -0.979
## potass 0.329 -0.942
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
## Cumulative Var 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000
summary(pc_S)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 41.2877486 3.880212624 2.64197339 1.624583979
## Proportion of Variance 0.9856182 0.008705172 0.00403574 0.001525986
## Cumulative Proportion 0.9856182 0.994323381 0.99835912 0.999885108
## Comp.5 Comp.6 Comp.7 Comp.8
## Standard deviation 0.353951757 2.561722e-01 8.510631e-02 2.372715e-02
## Proportion of Variance 0.000072436 3.794288e-05 4.187837e-06 3.255049e-07
## Cumulative Proportion 0.999957544 9.999955e-01 9.999997e-01 1.000000e+00
pc_R$loadings
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## rblood -0.194 0.417 0.400 0.652 0.175 -0.363 0.176 0.102
## plate -0.400 0.154 0.168 -0.848 0.230 -0.110
## wblood -0.459 0.168 -0.274 0.251 0.403 0.677
## neut -0.430 -0.472 -0.171 0.169 0.118 -0.237 0.678
## lymph 0.494 0.360 -0.180 -0.139 0.136 0.157 0.724
## bilir 0.319 -0.320 -0.277 0.633 -0.162 0.384 0.377
## sodium 0.177 -0.535 0.410 -0.163 -0.299 -0.513 0.367
## potass 0.171 -0.245 0.709 0.198 0.469 -0.376
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
## Cumulative Var 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000
summary(pc_R)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 1.6710100 1.2375848 1.1177138 0.88227419 0.78839505
## Proportion of Variance 0.3490343 0.1914520 0.1561605 0.09730097 0.07769584
## Cumulative Proportion 0.3490343 0.5404863 0.6966468 0.79394778 0.87164363
## Comp.6 Comp.7 Comp.8
## Standard deviation 0.69917350 0.66002394 0.31996216
## Proportion of Variance 0.06110545 0.05445395 0.01279697
## Cumulative Proportion 0.93274908 0.98720303 1.00000000
Las componentes principales, por construcción, están hechas para explicar el mayor porcentaje de la varianza. Por lo tanto, cuando se utiliza \(S\), como una variable (plate) tiene muchísima más variabilidad que las demás, la primera componente principal le asignará un coeficiente muy grande.
Aún más notorio es al comparar la proporción de la varianza que explica cada componente: cuando se utilizan las covarianzas, la primer componente principal explica casi toda. Peor: para explicar la misma proporción con el modelo de correlaciones se necesitan siete de las ocho componentes.
El análisis usando \(S\) sólo dice lo que podíamos ver desde \(\mathbf{\sigma}\): que plate varía mucho, y en este sentido es mejorel de \(R\). Sin embargo, hay que ser cuidadosos al interpretar porque los de \(R\) son de las variables estandarizadas, que están en desviaciones desde la media y no las medidas originales.
Considere los datos sobre venta de casas en T7-1.DAT. Este archivo tiene las siguientes variables:
Estime el modelo de regresión lineal \(\mathbf{Y} = \mathbf{\beta'x}+\mathbf{\varepsilon}\). Haga el diagnóstico del modelo sobre los residuales.
casas <- read_csv('./data/t71.csv')
casas_lm <-lm(p_venta~tamaño+p_valuacion, data=casas)
casas_tidy <- tidy(casas_lm)
casas_glance <- glance(casas_lm)
Parece ser que la valuación anterior no es estadísticamente significativa. Además, \(R^2=0.8344\) y \(R^2_{\textrm{aj}}=0.8149\), por lo que el modelo explica buena parte de la varianza.
Primero queremos buscar no-linealidad en los datos a través de una gráfica de residuales contra valores ajustados:
studentized <- rstudent(casas_lm)
casas_lm <- augment(casas_lm) %>%
mutate(.student.resid = studentized)
rm(studentized)
casas_lm %>%
ggplot(aes(x=.fitted, y=.resid)) +
geom_point() +
geom_smooth() +
geom_smooth() +
labs(x='Fitted values', y='Residuals')
casas_lm %>%
ggplot(aes(x=.fitted, y=sqrt(abs(.std.resid)))) +
geom_point() +
geom_smooth() +
ggtitle('Scale-location') +
theme(plot.title = element_text(hjust = 0.5)) +
labs(x='Fitted values', y=expression(sqrt(.std.resid)))
Aunque la LOESS parece tener forma de u en cerca de 75,pero en general la dispersión de los puntos no parece seguir patrones y la varianza se ve más o menos constante.
Queremos concentrarnos en un valor que puede ser problemático porque se ve lejos de los demás. Aprovechando que son sólo dos variables, podemos hacer una gráfica de dispersión para corroborar que efectivamente está lejos de los demás puntos.
casas_lm %>%
ggplot(aes(x=tamaño, y=p_valuacion)) +
geom_point()+
ggtitle('Predictor space') +
theme(plot.title = element_text(hjust = 0.5))
Para medir la influencia del punto aunque sea visualmente,
casas_lm %>%
ggplot(aes(x=1:20, y=.cooksd)) +
geom_point() +
geom_line() +
labs(x='Observation number', y="Cook's Distance")
casas_lm %>%
ggplot(aes(x=1:20, y=.hat)) +
geom_point() +
geom_line()+
labs(x='Observation number', y='Leverage')
p <- casas_lm %>%
ggplot(aes(x=.hat, y=.student.resid)) +
geom_point(aes(size=.cooksd)) +
labs(title='Influence plot',
x='Leverage',
y='Studentized residuals',
size="Cook's Distance") +
theme(plot.title=element_text(hjust=0.5))
ggplotly(p)
La distancia de Cook indica que la observación 16 podría tener mucha influencia sobre el modelo, y el apalancamiento confirma que es por estar lejos de las demás en el espacio de variables. Sin embargo, en la gráfica de influencia no se ve que el el valor de la respuesta esté muy lejos de la tendencia.
Para buscar correlación entre los errores,
casas_lm %>%
ggplot(aes(x=1:20, y=.resid)) +
geom_line() +
geom_point() +
labs(x='Observation number', y='Residuals')
Acf(casas_lm$.resid)
Y parece que no hay correlación significativa entre errores. Para verificar que se distribuyan normal y poder hacer inferencia tranquilamente, hacemos una qqplot
casas_lm %>%
ggplot(aes(sample=.student.resid)) +
geom_qq_band(alpha=0.5, distribution='norm') +
#geom_qq_band(alpha = 0.5, distribution='t', dparams=17) +
stat_qq_line() +
stat_qq_point() +
labs(x = "Theoretical quantiles", y = "studentized residuals")
Aunque los puntos no forman una línea muy recta en las colas, sólo uno se sale de los intervalos de confianza, lo cual es normal para una prueba al 5% de significancia sobre 20 datos. Nos damos por satisfechos suponiendo normalidad.
x0 <- c(1, 17, 46)
X <- casas %>%
mutate(const=1) %>%
select(const, tamaño, p_valuacion) %>%
as.matrix()
Y <- casas$p_venta
s2 <- casas_glance$sigma^2
beta_gorro <- casas_tidy$estimate
exts <- qt(0.95/2, 20-2-1, lower.tail = FALSE)*sqrt(s2*(1+t(x0)%*%solve(t(X)%*%X)%*%x0))
x0%*%beta_gorro-exts
## [,1]
## [1,] 77.41626
x0%*%beta_gorro+exts
## [,1]
## [1,] 78.24337
Realice una prueba LRT sobre $ H_0 : _2 = 0$ con nivel de significancia de \(\alpha= 0.05\)
modelo_H0 <- lm(p_venta ~ tamaño, data=casas)
modelo_Ha <- lm(p_venta ~ tamaño+p_valuacion, data=casas)
anova(modelo_H0, modelo_Ha)
## Analysis of Variance Table
##
## Model 1: p_venta ~ tamaño
## Model 2: p_venta ~ tamaño + p_valuacion
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 18 205.3
## 2 17 205.0 1 0.3027 0.0251 0.876
glance(modelo_Ha)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## * <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.834 0.815 3.47 42.8 2.30e-7 3 -51.7 111. 115.
## # ... with 2 more variables: deviance <dbl>, df.residual <int>
F_stat <- ((205.3-205)/(18-17))/glance(modelo_Ha)$sigma^2 #(Leí los numeritos de las llamadas de arriba)
p_value <- 1-pf(F_stat, 1, 17)
p_value
## [1] 0.8765281
No hay evidencia para rechazar la hipótesis nula, por lo que nos quedamos con el modelo que no incluye la valuación pasada.
Obtenga la gráfica de la \(C_p\) de Mallows para este problema.
cp <- regsubsets(p_venta~tamaño+p_valuacion, data=casas) %>%
summary() %>%
'$'(cp)
data_frame(p=1:2, cp=cp) %>%
ggplot(aes(p,cp)) +
geom_point() +
geom_abline(aes(slope=1,intercept=0))
Una vez más, indica que nos quedemos con el modelo más pequeño.
El Berkeley Guidance Study fue un estudio hecho a través del tiempo para seguir a un grupo de niños y niñas desde que nacieron en Berkeley entre enero de 1928 y junio de 1929, hasta al menos los 18 años. Los datos se encuentran en el archivo BGSall.DAT.
Para las niñas, obtengn las estadísticas sumarias usuales (medias, desviaciones estándar y correlaciones) para todas las variables, excepto Case y Sex. Obtenga la matriz de scatterplots para las variables de la edad 2, las variables de edad 18 y HT18. Resuma la información que se pueda extraer de la gráfica.
colnames <- c('Sex', 'WT2', 'HT2', 'WT9', 'HT9', 'LG9', 'ST9',
'WT18', 'HT18', 'LG18', 'ST18', 'soma', 'case')
berkeley <- read.table('./data/BGSall.DAT', skip=24, col.names=colnames) %>%
as_data_frame()
rm(colnames)
women_mean <- berkeley %>%
filter(Sex==1) %>%
select(-Sex, -soma) %>%
colMeans()
women_cor <- berkeley %>%
filter(Sex==1) %>%
select(-Sex, -soma) %>%
cor()
women_sd <- berkeley %>%
filter(Sex==1) %>%
select(-Sex, -soma) %>%
as.matrix() %>%
apply(2, sd)
ggpairs(data = berkeley,
mapping=ggplot2::aes(colour=as.factor(Sex)),
columns=2:12,
title='Berkeley Study Group'
)
Hay correlación positiva entre HT18 y casi todas las variables, particularmente la estatura a los dos años y a los nueve.
Ajuste el modelo de regresión lineal \[ \textrm{HT18}|\mathbf{X} = \alpha_0 + \alpha_1\textrm{WT2} + \alpha_2\textrm{HT2} + \alpha_3\textrm{WT9} + \alpha_4\textrm{HT9} + \alpha_5\textrm{LG9} + \alpha_6\textrm{ST9} + \epsilon \] suponiendo que la varianza es constante y dé los estimados y los errores estándar para todos los parámetros, así como el valor del coeficiente de determinación.
bk_model <- lm(HT18~WT2+HT2+WT9+HT9+LG9+ST9, data=berkeley)
tidy(bk_model)
## # A tibble: 7 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 41.3 21.7 1.91 0.0590
## 2 WT2 1.42 0.494 2.87 0.00478
## 3 HT2 0.127 0.255 0.496 0.621
## 4 WT9 -0.0898 0.270 -0.333 0.740
## 5 HT9 0.969 0.177 5.46 0.000000236
## 6 LG9 -1.20 0.569 -2.12 0.0360
## 7 ST9 0.0994 0.0433 2.29 0.0235
glance(bk_model)$adj.r.squared
## [1] 0.5124361
Muestre numéricamente que \(R^2\) es el mismo valor que el coeficiente de correlación entre HT18 y los valores ajustados de la ecuación de la pregunta anterior. Haga una gráfica de la respuesta contra los valores ajustados, y dé una interpretación visual del coeficiente de determinación.
r_squared <- glance(bk_model)$r.squared
bk_model_aug <- bk_model %>%
augment()
corr_calc <- cor(bk_model_aug$HT18, bk_model_aug$.fitted)^2
r_squared/corr_calc
## [1] 1
bk_model_aug %>%
ggplot(aes(HT18, .fitted)) +
geom_point() +
geom_smooth(method='lm', se = F)
El valor de \(R^2\) dice que la recta explica cerca del 53% de la variabilidad en las observaciones, cosa que vemos en que los puntos están muy dispersos al rededor de la línea.
De las unidades de medición de cada uno de los coeficientes de regresión estimados.
Obtenga pruebas de que cada coeficiente es igual a 0, y dar el valor de la estadística de prueba, su p−valor, y un breve resumen del resultado.
Haciendo las pruebas margianles:
summary(bk_model)
##
## Call:
## lm(formula = HT18 ~ WT2 + HT2 + WT9 + HT9 + LG9 + ST9, data = berkeley)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.4766 -4.5606 0.1421 4.0278 11.7021
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.29577 21.67489 1.905 0.05898 .
## WT2 1.41695 0.49351 2.871 0.00478 **
## HT2 0.12656 0.25501 0.496 0.62054
## WT9 -0.08979 0.26981 -0.333 0.73984
## HT9 0.96851 0.17740 5.460 2.36e-07 ***
## LG9 -1.20494 0.56873 -2.119 0.03604 *
## ST9 0.09936 0.04333 2.293 0.02346 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.176 on 129 degrees of freedom
## Multiple R-squared: 0.5341, Adjusted R-squared: 0.5124
## F-statistic: 24.65 on 6 and 129 DF, p-value: < 2.2e-16
Los valores p indican que sólo podemos rechazar la hipótesis de coeficientes nulos con confianza superior a 95% en WT2, HT9, LG9 y ST9.
Probar la hipótesis de que los predictores a edad 2 no son necesarios en el modelo (i.e. prueba que \(NH : \alpha_1 = \alpha_2 = 0\) versus la alternativa general).
anova(lm(HT18~WT9+HT9+LG9+ST9, data=berkeley),
bk_model)
## Analysis of Variance Table
##
## Model 1: HT18 ~ WT9 + HT9 + LG9 + ST9
## Model 2: HT18 ~ WT2 + HT2 + WT9 + HT9 + LG9 + ST9
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 131 5380.8
## 2 129 4919.8 2 460.99 6.0437 0.003098 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La estadística F da suficiente evidencia para rechazar la hipótesis nula con 97% de confianza, indicando que es mejor el modelo con las variables a los dos años.
De intervalos de confianza del 95 % para \(\alpha_4\) y para \(\alpha_4 - \alpha_2\) de la ecuación.
Suponiendo normalidad y homocedasticidad, los coeficientes siguen \[ \mathbf{\hat{\alpha}} \sim \mathcal{N}(\mathbf{\alpha}, \sigma^2(X^\top X)^{-1}) \]
Por lo que podemos dar intervalos simultáneos
alpha <- bk_model$coefficients
s_alpha <- vcov(bk_model)
print('Intervalo para alpha 4')
## [1] "Intervalo para alpha 4"
lim_int <- sqrt(s_alpha[5,5])*sqrt((6+1)*pf(6+1,136-(6+1),0.05, lower.tail=F))
alpha[[4]] + c(-1,1)*lim_int
## [1] -0.5302151 0.3506350
print('')
## [1] ""
print('Intervalo para alpha 4 - alpha 2')
## [1] "Intervalo para alpha 4 - alpha 2"
a <- c(0, 0, -1, 0, 1, 0, 0)
centro <- a%*%alpha
lim_int <- sqrt(t(a)%*%s_alpha%*%a)* sqrt(((6+1)*(136-1)/(136-(6+1)))*pf(6+1,136-(6+1),0.05, lower.tail=F))
c(centro -lim_int, centro+lim_int)
## [1] -0.1317404 1.8156470
Con datos de alguna fuente de datos abiertos, proponer (a) un modelo de regresión múltiple, y (b) un modelo de componentes principales. En cada caso, realizar el análisis correspondiente, y elaborar un reporte del análisis en cada caso, de no más de una cuartilla en cada caso.
rm(list=ls())
power_plant <- read_excel('./data/CCPP/Folds5x2_pp.xlsx')
names(power_plant) <- str_replace_all(names(power_plant), ' ', '_')
pp_models <- regsubsets(PE~., data=power_plant, nvmax = 2500)
data_frame(p=1:4, cp=summary(pp_models)$cp) %>%
ggplot(aes(p,cp)) +
geom_point() +
geom_line()
# Decidimos quedarnos con el de tres
summary(pp_models)
## Subset selection object
## Call: regsubsets.formula(PE ~ ., data = power_plant, nvmax = 2500)
## 4 Variables (and intercept)
## Forced in Forced out
## AT FALSE FALSE
## V FALSE FALSE
## AP FALSE FALSE
## RH FALSE FALSE
## 1 subsets of each size up to 4
## Selection Algorithm: exhaustive
## AT V AP RH
## 1 ( 1 ) "*" " " " " " "
## 2 ( 1 ) "*" " " " " "*"
## 3 ( 1 ) "*" "*" " " "*"
## 4 ( 1 ) "*" "*" "*" "*"
pp_model_cp <- lm(PE~AT+V+RH, data=power_plant)
anova(lm(PE~., data=power_plant))
## Analysis of Variance Table
##
## Response: PE
## Df Sum Sq Mean Sq F value Pr(>F)
## AT 1 2505095 2505095 120563.32 < 2.2e-16 ***
## V 1 46766 46766 2250.72 < 2.2e-16 ***
## AP 1 6259 6259 301.23 < 2.2e-16 ***
## RH 1 29875 29875 1437.81 < 2.2e-16 ***
## Residuals 9563 198702 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pp_model_anova <- lm(PE~., data=power_plant)
glance(pp_model_cp)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC
## * <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 0.928 0.928 4.57 41321. 0 4 -28110. 56229.
## # ... with 3 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>
glance(pp_model_anova)
## # A tibble: 1 x 11
## r.squared adj.r.squared sigma statistic p.value df logLik AIC
## * <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 0.929 0.929 4.56 31138. 0 5 -28088. 56188.
## # ... with 3 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>
anova(pp_model_cp, pp_model_anova)
## Analysis of Variance Table
##
## Model 1: PE ~ AT + V + RH
## Model 2: PE ~ AT + V + AP + RH
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 9564 199598
## 2 9563 198702 1 895.28 43.087 5.507e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ambos modelos tienen \(R^2\) y AIC similares, pero la prueba F declara mejor al modelo completo.
pp_model_data <- pp_model_anova %>%
augment()
pp_model_data %>%
ggplot(aes(x=.fitted, y=.resid)) +
geom_point(alpha=0.5) +
geom_smooth()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Se alcanza a distinguir que en valores pequeños los residuales son mayores.
pp_model_data %>%
ggplot(aes(1:9568, .cooksd)) +
geom_point() +
geom_line()
pp_model_data %>%
ggplot(aes(1:9568, .hat)) +
geom_point() +
geom_line()
Ningún punto parece tener demasiado alto apalancamiento.
pp_model_data %>%
ggplot(aes(sample=.std.resid)) +
geom_qq_band(alpha=0.5, distribution='norm') +
#geom_qq_band(alpha = 0.5, distribution='t', dparams=17) +
stat_qq_line() +
stat_qq_point() +
labs(x = "Theoretical quantiles", y = "standarized residuals")
Pero los errores definitivamente no son normales, por lo que no podemos hacer inferencia sobre los coeficientes.
Haciendo la descomposición en componentes principales
pp_pc <- power_plant %>%
select(-PE) %>%
prcomp(scale=T) %>%
print()
## Standard deviations (1, .., p=4):
## [1] 1.5615895 0.9533469 0.7416225 0.3202560
##
## Rotation (n x k) = (4 x 4):
## PC1 PC2 PC3 PC4
## AT 0.6148135 -0.04980763 0.1667322 -0.7692359
## V 0.5596829 0.09961674 0.5936655 0.5695549
## AP -0.4041859 -0.64834308 0.6287465 -0.1447856
## RH -0.3813045 0.75315799 0.4737504 -0.2508395
summary(pp_pc)
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.5616 0.9533 0.7416 0.32026
## Proportion of Variance 0.6096 0.2272 0.1375 0.02564
## Cumulative Proportion 0.6096 0.8369 0.9744 1.00000
Aunque la primera componente principal es difícil de interpretar, la segunda tiene una interpretación que coincide con lo dicho en la documentación de los datos: separa vacumm —que afecta a las turbinas de vapor— , de las tres variables ambientales, que afectan a las de gas.
Aunque se necesitan tres de las cuatro para explicar la variabilidad de los predictores, vamos a ver cuántas se necesitan para explicar la de la respuesta.
pcr_model <- pcr(PE~., data=power_plant)
summary(pcr_model)
## Data: X dimension: 9568 4
## Y dimension: 9568 1
## Fit method: svdpc
## Number of components considered: 4
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps
## X 62.99 92.27 98.41 100.00
## PE 62.17 84.83 85.79 92.87
Con la variabilidad de la respuesta nos va peor.
predplot(pcr_model)
validationplot(pcr_model, val.type = 'MSEP')
Sin embargo, en términos del error medio es parece que quedarnos con dos componentes es lo mejor. Hacemos la prueba formal
predictores_cp <- pp_pc$x %>%
as_data_frame() %>%
mutate(PE=power_plant$PE) %>%
select(PE, everything()) %>%
print()
## # A tibble: 9,568 x 5
## PE PC1 PC2 PC3 PC4
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 463. -1.67 -1.25 0.449 -0.339
## 2 444. 0.747 -1.44 0.784 -0.104
## 3 489. -2.27 1.07 -0.527 0.536
## 4 446. 0.351 0.517 -0.0437 0.0267
## 5 474. -1.80 1.57 -0.653 -0.144
## 6 444. 1.22 -0.642 -0.193 -0.178
## 7 467. -0.868 -0.0394 -0.424 -0.127
## 8 478. -1.48 -1.00 -0.279 0.595
## 9 476. -0.566 -2.62 -0.685 0.443
## 10 478. -1.19 -0.370 -0.564 0.334
## # ... with 9,558 more rows
names(predictores_cp) <- c('PE', 'PC1', 'PC2', 'PC3', 'PC4')
pcr_model_full <- lm(PE~.,data=predictores_cp)
summary(pcr_model_full)
##
## Call:
## lm(formula = PE ~ ., data = predictores_cp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.435 -3.166 -0.118 3.201 17.778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 454.36501 0.04660 9750.14 <2e-16 ***
## PC1 -9.99355 0.02984 -334.87 <2e-16 ***
## PC2 -1.53914 0.04888 -31.49 <2e-16 ***
## PC3 -5.08334 0.06284 -80.89 <2e-16 ***
## PC4 10.16892 0.14552 69.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.558 on 9563 degrees of freedom
## Multiple R-squared: 0.9287, Adjusted R-squared: 0.9287
## F-statistic: 3.114e+04 on 4 and 9563 DF, p-value: < 2.2e-16
pcr_model_full %>%
augment() %>%
ggplot(aes(sample=.std.resid)) +
geom_qq_band(alpha=0.5, distribution='norm') +
stat_qq_line() +
stat_qq_point() +
labs(x = "Theoretical quantiles", y = "standarized residuals")
Lamentablemente, no podemos hacer inferencia sobre los coeficientes porque el error no es normal.